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Abstract. 

Much of standard galaxy dynamics rests on the implicit as- 
sumption that the corresponding iV-body problem is (near) in- 
tegrable. This notion although leading to great simplification is 
by no means a fact. In particular, this assumption is unlikely to 
be satisfied for systems which display chaotic behaviour which 
manifests itself on short time-scales and for most initial con- 
ditions. It is therefore important to develop and test methods 
that can characterize this kind of behaviour in realistic situa- 
tions. We examine here a method, pioneered by Krylov (1950) 
and first introduced to gravitational systems by Gurzadyan 
& Savvidy (1984,1986). It involves a metric on the configura- 
tion manifold which is then used to find local quantification of 
the divergence of trajectories and therefore appears to be suit- 
able for short time characterization of chaotic behaviour. We 
present results of high precision iV-body simulations of the dy- 
namics of systems of 231 point particles over a few dynamical 
times. The Ricci (or mean) curvature is calculated along the 
trajectories. Once fluctuations due to close encounters are re- 
moved this quantity is found to be almost always negative and 
therefore all systems studied display local instability to ran- 
dom perturbations along their trajectories. However it is found 
that when significant softening is present the Ricci curvature is 
no longer negative. This suggests that smoothing significantly 
changes the structure of the 6N phase space of gravitational 
systems and casts doubts on the continuity of the transition 
from the large- TV limit to the continuum limit. From the value 
of the negative curvature, evolution time-scales of systems dis- 
playing clear instabilities (for example collective instabilities 
or violent relaxation) are derived. We compare the predictions 
obtained from these calculations with the time-scales of the 
observed spatial evolution of the different systems and deduce 
that this is fairly well described. In all cases the results based on 
calculations of the scalar curvature qualitatively agree. These 
results suggest that future applications of these methods to re- 
alistic systems may be useful in characterizing their stability 
properties. One has to be careful however in relating the time- 
scales obtained to the time-scales of energy relaxation since 
different dynamical quantities may relax at different rates. 
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1. Introduction and motivation 

1.1. On the assumptions of standard galaxy dynamics 

The classical assumption often made in galaxy dynamics (Bin- 
ney & Tremaine 1987 (BT) chapter 4) is that present day 
galaxies can be treated as collisionless fluids in steady states 
described by one particle distribution functions obeying the 
time independent collisionless Boltzmann equation (CBE). In 
addition, it is often supposed that these distribution functions 
are completely characterized by the isolating integrals of mo- 
tion (of single stars). 

implies that almost all orbits are regular and therefore stars 
conserve as many integrals of motion (in involution) as the 
number of spatial dimensions they move in. Thus, for a three 
dimensional system of N stars there are 3N conserved quan- 
tities and the iV-body problem is solvable by quadratures — or 
integrable (see e.g., Whittaker 1937 or Goldstein 1980 for dis- 
cussions using classical analysis: more modern treatments can 
be found in Arnold 1989 (ARN) or Abraham & Marsden 1978). 

Thus in this picture galaxies are modeled as classical fluids 
with time independent densities which give rise to potentials 
similar to those for which the Hamilton- Jacobi equation is sep- 
arable. Obviously these assumptions cannot be strictly satisfied 
for real galaxies. It is argued however, that due to the large two 
body relaxation time (e.g., Saslaw 1985) 

6 32irG 2 m 2 n\n{N/2) 1 ' 

(where v is the characteristic speed of a star and n is the num- 
ber density of "field stars" of mass m) that these systems can 
be treated as effectively collisionless over many Hubble times. 
However all that Eq. (1) is saying is that the direct impulse 
from encounters between pairs of stars is small relative to the 
velocity of a star in the smoothed potential. This is due to the 
long range nature of gravitational interactions which ensures 
that the whole system contributes to the mean force on a star 
at all times. In the approximation leading to (1) however the 
perturbation to the trajectory of a test star due to a certain 
field star is only calculated at their closest approach — that 
is only once during a crossing time. Nevertheless, due to the 
complicated nature of the solutions of the Newtonian equations 
for N > 2, discreteness can have a major indirect effect; the 
nonlinearity of the problem prevents the adding up of the mo- 
tions due to individual binary interactions to each other and to 
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the motion in the mean field since although the forces add up 
linearly the solutions of the Newtonian equations do not. This 
suggests that the TV-body problem is perhaps best studied in 
its entirety. In this regard it may be perhaps useful to recall 
that the collisionless steady-state approximation is not a trivial 
simplification of the dynamics: It reduces the TV-body problem 
to TV one-particle problems in a given potential, thus reducing 
the number of coupled first order differential equations to be 
solved from 6TV equations to 6 equations. 

For general Hamiltonian dynamical systems, regular mo- 
tions form invariant (under time propagation of the solution) 
tori in the 6TV phase-space. In integrable systems these oc- 
cupy the whole of that space. Under perturbations, however, 
some tori are destroyed leaving behind volumes of phase space 
where irregular (or chaotic) motion can occur. The extent of 
the chaotic region will depend on the strength of the perturba- 
tion. Bounds given by the KAM theorem predict that a posi- 
tive measure of tori will survive provided the perturbation to 
the potential /i < fi c , where /j, c is a critical amplitude (e.g., 
Mackay & Meiss 1987). In general, the value of /i c decreases as 
exp(— TV log TV) and has been shown to be irrelevant for many 
higher dimensional physical systems (Pettini 1993 (P93) and 
the references therein) while the perturbation due to discrete- 
ness noise decreases only as l/y/N if it is random (e.g., Saslaw 
1985). The diffusion time away from the remaining tori is also 
a decreasing function of TV; td ~ exp(l//i) 1,/JV (Nekhoroshev 
1977; Perry & Wiggins 1994). These results suggest that a sys- 
tem is unlikely to become more regular with increasing TV, in 
direct contradiction with the results of two body relaxation 
estimates but as expected from considering the complicated 
nature of solutions of generic TV-body systems. In fact, for 
the particular case of a spherical gravitational system, it was 
shown by Gurzadyan & Savvidy (1986) (GS) that as TV in- 
creases this system tends towards an Anosov (1967) C-system 
with maximal instability in phase-space (if one neglects escapes 
and direct collisions) . Making the same kind of assumptions as 
those used to obtain Eq. (1) (by considering an infinite and 
homogeneous medium), GS obtain the following relation for 
the TV-body relaxation time arising from the instability 
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which is considerably shorter than the binary relaxation time. 

C-systems are very irregular — to the point that their evo- 
lution can be described as a Markow process (e.g., Pesin 1989). 
The collisionless approximation however predicts that spherical 
systems end up in completely integrable steady states — the 
most regular and predictable systems that can exist. The con- 
tradiction can possibly arise because the analysis leading to (1) 
focuses on the fact that the force function becomes smoother 
as TV increases. It then assumes directly that this implies the 
solution becoming more regular not taking into account that 
exponentially smaller perturbations arc sufficient with increas- 
ing TV to make an TV-body system unstable. Therefore, while 
some quantities that do not depend on the exact details of the 
dynamics may be slowly evolving and not sensitively dependent 
on the discreteness noise, the solutions themselves are heavily 
sensitive to noise. For example, the change in the energies of 
stars that arises from discreteness noise is likely to be much 
slower than the changes in their trajectories. A rough exam- 
ple is the fact that numerically integrated orbits in fixed (and 
smooth) potentials can have the energy (and all the Poincare 



invariants) conserved along their trajectories to up to ten dig- 
its while the trajectories themselves are completely inaccurate 
perhaps not resembling (even qualitatively) the real ones (El- 
Zant 1996b). The numerical errors here represent the noise (if 
they are taken to be random, a proposition still under discus- 
sion: McCauley 1993). This is because here energy is a scalar 
function of the phase space variables and is related to the force 
by a path independent integral. That path on the other hand 
is a 6 component vector in phase space. One can thus easilly 
envisage perturbations that change the trajectory of a parti- 
cle without changing its energy. A long time-scale of energy 
relaxation does not therefore imply that the actual detailed 
dynamics are collisionless in the sense of being stable solutions 
of the CBE. One would expect that some macroscopic quan- 
tities could be affected, for example velocity dispersion or the 
actual shape of the gravitational system. 

In fact, simulations of single particle motions in fixed po- 
tentials show that the time-averaged phase-space density dis- 
tribution of individual chaotic trajectories, as well as their sta- 
tistical properties, change significantly over a time-scale <C r b 
when the potential is given some graininess (Pfenniger 1986; 
Udry & Pfenniger 1988), or when the trajectories are given 
small random kicks without changing their total energies signif- 
icantly (Kandrup 1994; Merritt & Valluri 1996; El-Zant 1996b). 
The full TV-body problem would be expected to be much more 
irregular and prone to evolution on small (compared to r 6 ) 
time-scales. Such effects are actually observed in TV-body sim- 
ulations of up to 10 6 heavily softened particles but are assumed 
to be due to the relatively small number of particles present, 
even though the two body relaxation time in these simulations 
is still usually much larger than the time-scales considered, (a 
review of such occurrences is given by Hernquist & Ostriker 
1992: see also van Albada 1986; Sellwood 1987: Zhang 1996 
presents a detailed study of an example of a case where dis- 
creteness noise interacts with global irregularities in the density 
thus triggering evolution). In addition, processes involving non- 
stellar objects such as interactions with giant molecular clouds 
in galactic disks or black holes in halos, small dissipative per- 
turbations (e.g., Pfenniger & Norman 1990) etc. can cause even 
more serious trouble for the collisionless approximation. 

Even if the collisionless approximation does hold, this does 
not guarantee that a given density distribution would not 
evolve over a Hubble time. This is because individual trajec- 
tories, even in a smooth steady-state potential, can have time 
dependent density distributions over such a time-scale. This 
would be the case in general non-spherical potentials as noted 
by Binney (1982). Hasan et al. (1993) describe such behaviour 
for barred spirals while Merritt & Fridman (1996) show that 
it may also be important for ellipticals. El-Zant (1996b) exam- 
ines the case of trajectories started near the symmetry plane 
of disks embedded in triaxial halos. Finally, it seems that per- 
haps there is also some observational evidence for evolutionary 
phenomena occuring in galaxies (Wielen 1977; Pfenniger et al. 
1994; Courteau et al. 1996; Pucacco 1992). 

The above considerations suggest (although in no way 
prove) that even though gravitational systems obey the CBE 
in the infinite TV limit (e.g., Braun & Hepp 1977) the exis- 
tence and stability (against discreteness noise) of steady state 
solutions of that equation are in question when a significant 
amount of chaos is present. This means that while in some 
cases the classical theory may still hold over a few Gyr (con- 
sidering the primitive evolutionary state of cold disks in some 
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spiral galaxies one has to admit that this could often be the 
case) it does not always obviously do so. 

1.2. How chaos drives evolution 

A distribution of stars giving rise to an integrable potential 
can either oscillate coherently or reach a macroscopic steady 
state through phase mixing. This is a trivial form of relax- 
ation which is simply due to stars moving at different angular 
frequencies on their respective KAM tori in the 6-dimensional 
phase-space — the motion in the 6iV-dimensional phase-space 
being a 3iV torus characterized by the 3JV integrals. Phase 
mixing conserves the action variables characterising this torus 
and which are crucial in determining the physical state of a sys- 
tem. Most galaxies are not believed to be undergoing significant 
large scale oscillations, we therefore conclude that if they are 
integrable dynamical systems they must be in a steady state. 
Moreover, if they are sufficiently far from any chaotic system 
they must be stable to small perturbations (KAM theory). This 
is what is assumed in the classical theory. 

A necessary condition for evolution therefore is the non- 
integrability of the system. This condition is satisfied for N- 
body gravitational systems since there are no global integrals of 
motion other than the classical ones (Poincare 1889). However, 
not all non-integrable systems exhibit interesting behaviour 
(different from classical theory) in the time-scales of interest. 
We need the property of phase-space mixing — that is the 
spread of localized volume elements (corresponding to sets of 
initial conditions) to cover large areas of the 6N phase-space 
(while still conserving their original Lebesgue measure: e.g., 
Sagdeev et al. 1988). Phase-space mixing leads to diffusion in 
the action variables and may therefore lead to evolution in a 
system's physical parameters. For this process to be efficient 
it is necessary that the system be sufficiently chaotic so that 
the diffusion occurs over short enough time-scales and covers 
a large range of initial states, only then can we say that the 
macroscopic state corresponding to a given set of micro-states 
can evolve over a Hubble time. 

It is therefore important to examine different methods 
for detecting these processes and the time-scales associated 
with them in practical situations where the predictions can be 
tested. The rest of this study is a step in this direction. In the 
next section we describe why certain subtleties related to N- 
body gravitational systems require local methods to be used 
in the quantification of chaotic behaviour. One such approach 
based on the geometry of the configuration manifold is then 
described. In section 3 we describe some of the possible ap- 
plications of this approach while in section 4 results of some 
numerical experiments testing the method are reported. 

2. Characterization of chaos and the Ricci criterion 

2.1. Difficulties with commonly used methods 

Central to the idea of phase-space mixing is the concept of 
dynamical entropy used to quantify it. The most important of 
such quantities is the so called Kolmogorov-Sinai entropy. The 
easiest way of calculating the KS entropy is by evaluating the 
Liapunov exponents, which for a dynamical system defined by 
the vector equations 



^o,X ) = £milogj[||l, 
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where £ is a tangent space vector arising from the solution of 
the variational equations 

i = D x F(X(M ,X ))£. (5) 
The KS entropy is then usually given by (e.g., LL) 



KS ■ 



a(e,X )dX , 
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where the sum is taken over all positive exponents and the in- 
tegral is over all possible initial conditions. For systems with 
simple enough phase-space and when the infinite time charac- 
teristics are required, the KS entropy calculated with the help 
of this formula suffices to describe the ergodic properties of 
a system. Positive KS entropy over a compact phase-space, or 
some region of it that has this property, is a sufficient condition 
for the presence of what is usually called chaos and the accom- 
panying erratic behaviour which leads to the approach towards 
statistical equilibrium even in low dimensional systems (e.g., 
McCauley 1993). The evolution time-scale is usually related to 
the inverse of the KS entropy. 

Open systems interacting via un-softened Newtonian po- 
tentials do not have compact phase-spaces however, and there- 
fore the situation is more complicated. Here, no final state 
exists and one has to distinguish between the various stages 
of evolution i) Violent relaxation, ii) Collective (plasma type) 
instabilities iii) Evolution towards an isotropic rotator and iv) 
Kinetic evolution towards equipartition of energy, core collapse 
etc. Here, the distinctions are rather arbitrary and are used 
only for clarity — these processes may interact and influence 
each other. Stage three has been given very little attention in 
the literature because of the lack of a mechanism from tradi- 
tional physics leading to such evolution (on a time-scale smaller 
than given by (1)). Nevertheless, the chaotic nature of the N- 
body problem may provide a clue as to how this might happen. 

All the above processes should be characterized by phase- 
space instability leading to mixing, but evidently cannot be de- 
scribed by any quantities defined only for infinite times. One 
way out of this problem is to integrate either the linearized 
equations (5) (e.g., Goodman et al. 1993) or the full non-linear 
equations (3) for slightly different initial conditions (e.g., Kan- 
drup et al. 1994 and the references therein) for short times. 
One drawback of such an approach, however, is that one is 
comparing the divergence in phase-space of different tempo- 
ral states and not trajectories. Integrable systems usually have 
linear (in time) phase-space divergence between neighbouring 
states — but only on average. A simple illustration of the type 
of problems involved is provided by examining the behaviour 
of a pendulum 



with linearized equation 
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X = F(X), 

with solution X = X(t, t , X Q ) are given by 



(3) When cos 9 is negative this has a solution 
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That is, during this interval, the solution of the linearized equa- 
tion predicts "exponential instability". Obviously, in this ex- 
ample, the trajectories do not diverge at all — but the states 
did. In addition, methods that evaluate the whole set of Lia- 
punov exponents of a dynamical system are fairly sophisticated 
(e.g., Eckmann & Ruelle 1985) and are not practical for higher 
dimensional systems. This in practice will mean that one will 
have to evaluate only the largest exponent which is a measure 
of the maximal instability in phase-space. Now, in a multidi- 
mensional separable system it is likely that at any given time 
there will be some oscillations that are in the "cos 6 < region" 
— that is displaying exponential divergence in the linearized 
dynamics. An initial randomly oriented vector in the linear 
tangent space will always reorient itself along the direction of 
maximum expansion under the stretching effect of the flow 
(e.g., Wolf et al. 1985). It could therefore be possible to obtain 
average exponential divergence in the linearized dynamics even 
for separable systems. This approach therefore may be unable 
to distinguish between genuine phase space mixing leading to 
evolution and trivial phase mixing of temporal states. It may 
also be useful to note here that there are few strict results con- 
cerning the properties of the Liapunov exponents for higher 
dimensional systems (Eckmann & Ruelle 1985; Pesin 1989). 

2.2. The study of motion on Lagrangian manifolds 

There are a variety of ways of transforming Hamiltonian prob- 
lems into the study of some metric space (see P93; Gurzadyan 
& Kocharyan 1994 or the articles by Gurzadyan and Pettini in 
Gurzadyan & Pfenniger 1994). The oldest and most well known 
of these hinges on the observation, apparently first made by 
Hertz (1900) in the course of his remarkable reformulation of 
classical mechanics, that the Maupertuis principle 



ds 



S / 2Tdt = 0, 



(10) 



(where T is the kinetic energy along the motion on a trajectory 
7) from which the equations of motion in their Lagrangian form 
arise is actually an expression for geodesies on a the configura- 
tion manifold M where the motion is restricted as a result of 
conservation laws. The fact that (10) defines geodesies is clear 
from the following relations: 

S e J Tdt = 8 e j VTVf dt = 8 e J VTdl = 8 e J VE-Vdl,{ll) 
where E = T + V is the total energy and dl = \J ciijdqv with 



dq i dq j 
' dt dt 



(12) 



for particles of unit mass. Here the elements of the 

metric tensor and 8 e refers to variations in the trajectory 7 
holding the energy and the end points fixed. These are thus 
geodesies in the energy sub- manifold of the configuration space 
and the q % are coordinates on it. If we now choose Cartesian 
coordinates in the enveloping 3N space then 



di 2 = J2( dxa ) 2 - 

The metric then is 



(13) 



(B-V)^(da 



(14) 



From the Jacobi equation which describes the geodesic devia- 
tion on M (e.g., Misner et al. 1973) 



V u V u n + Riem(n, u)u = 0, 



(15) 



(where n is a separation vector analogous to £ in (4), Riem is 
the Riemann curvature operator, and V u the covariant deriva- 
tive) one can obtain the following relation for the norm of the 
normal component of this deviation 



d 2 



n 



ds 2 
where 



+ 211 V u n 



[Riem(n, u)u] .n 



(16) 



(17) 



is the two dimensional curvature in a plane defined by u x n 
and where we have used the fact that nx.u = . If fc Ujn is 
negative everywhere for all planes as defined above (that is for 
all n normal to u) and if — k — min I &n,u 

I we have 



n(s) ||>i||n||exp< 



for n > 0, and 



n(s) ||< ± || n 



exp - 



(18) 



(19) 



for n < 0. These relations describe the linearized dynamics in 
the "dilating" and "contracting" spaces characteristic of the 
class of Anosov (1967) C-systems to which, as was mentioned 
earlier, large spherical JV-body systems belong. In this case, 
\/ — fc u ,n averaged over M is the KS entropy. In compar- 
ing two C-systems therefore one can define the system with 
larger average value of — ^2 fc Ujn to be more unstable. This 
system will have a larger exponentiation rate and so initial 
conditions will tend to mix faster along geodesies. To deter- 
mine how fast the initial conditions mix in time, we note that 
ds/dt = \p2T so that if T does not vary too much during the 
evolution s ~ \f2Tt. It is clear that a system with the above 
characteristics cannot conserve its action variables since there 
is always a deviation of trajectories normal to the motion in 
phase-space (which is the cotangent bundle of M). 

2.3. Ricci curvature and the corresponding criterion 

GS have shown that the condition for C-systems is not sat- 
isfied for general gravitational ones. However, as Kandrup 
(1990a,1990b) has shown, the probability of a two dimensional 
curvature being positive along a iV-body system's trajectory 
decreases exponentially with increasing N. Also, in the N- 
body problem relation (18) implies that all orbits are unsta- 
ble at all times. However, one needs much less than what is 
described by (18) for observable effects of instability to be de- 
tected (just 10% of orbit becoming unstable may be enough). 
On the other hand, just a few orbits being unstable in gen- 
eral would not significantly change the physical properties of 
a system. We therefore need some averaged form of (16) and a 
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corresponding instability relation instead of (18) to character- 
ize such behaviour. Ideally such a relation should not require 
the evaluation of the Riemann tensor. 

A natural way of proceeding is by using the Ricci (or mean) 
curvature of the manifold M 



r u (s) = Ri 



(20) 



where Rij are elements of the Ricci tensor and u % 



4r- are the 

as 



components of the geodesic velocity vector u. The Ricci curva- 
ture is related to the two dimensional curvatures by (Eisenhart 
1926) 



7"u — y ' fcn M ,u(s)- 



(21) 



The value of r u does not depend on the particular set of normal 
directions n chosen so that r u /(3iV — 1) can be seen as the 
average value of fc Ujn over all possible directions normal to u 
on the configuration manifold M. 

In the case when all k's are negative, y/—r u averaged over 
the whole manifold corresponds to the Kolmogorov entropy. In 
general, as was first noted by Gurzadyan & Kocharyan (1987), 
it will provide an "averaged" measure of irregularity. In these 
terms Equation (16) can be written as 



d 2 Z 2 
ds 2 



r-u(s) 



(3JV - 1) 



Z 2 + 2( 



V u n 



(22) 



where Z is now to be interpreted as the norm of a vector that is 
a member of a random field of vectors with uniform distribution 
in directions normal to u and equal magnitude (El-Zant 1996a). 

If r u (s) is negative on a region of M then one can obtain 
relations for Z analogous to those obtained for || n || 2 in (18) 
with —k r = ml "j^"i S ^ ■ One can then apply the criterion of rela- 
tive instability of C-systems to general Hamiltonian dynamical 
systems in regions of their configuration manifolds where the 
Ricci curvatures are negative, which in this case will express the 
relative probability of any two systems being unstable under 
random perturbations. We adopt here the following definition, 
convenient for numerical studies of iV-body systems. 

Definition: 

Let Ri be some subset of the configuration manifold Mi and R2 
be a subset of a manifold M2 with M2 not necessarily different 
from Mi. Suppose also that the Ricci curvature is negative 
in both of these regions. We will say that Ri corresponds to 
configurations of a dynamical system that are more unstable 
than those represented by R 2 if the average value of \/—r u is 
larger in Ri than in R 2 . 

Note that: As in the case of C-systems we obtain time- 
scales from the relation 4f = \/2T. If we are comparing systems 
with different kinetic energies, the evolutionary times derived 
will have to be scaled accordingly. If the kinetic energies of 
systems are changing in the region of the dynamical systems 
of interest then time-scales derived from the Ricci curvature 
alone are not rigorous. If the logarithmic time derivative of 
the kinetic energy is small however then s ~ Tt (where the 
bar denotes an average over the region of interest). Otherwise 
a fully dynamical formulation with time replacing s in Equa- 
tion (22) would have to be considered (P93; Cerruti-Sola & 
Pettini 1995). The latter approach would also have to be used 



if the curvature on M is not predominantly negative. Also if 
the systems we are comparing consist of different numbers of 
particles, r u will have to be divided by 3iV — 1 . 

The definition leaves us the choice to compare different 
areas of the manifold of the same dynamical system, or those 
of different systems. Also the averages can either be static, or 
taken along a computed trajectory of the system. Since the 
regions R in the definition above can be taken as small as we 
wish, the method is clearly local. This is possible because r u is 
directly related to the local geometry of the manifold where a 
system lives and is not an asymptotic quantity. Also, although 
the above formulation does not allow explicitly for dissipative 
forces, these can be added as time dependent perturbation to 
an open Hamiltonian system. 

To actually calculate the value of r u , one contracts the 
Riemann tensor (which for the metric (14) is given in GS) and 
uses (20) to obtain 



r„ = 3A- 



(W iU \ 
W 2 



~ 2A W 



W 3 



V 2 W 
' 2W 2 



(23) 



with A = 2^2. Here W denotes T = T(V) while VW and 
\7 2 W represent its Cartesian gradient and Laplacian respec- 
tively. From the metric (14) one can deduce that u 1 — d ^J^ 
and || u ||= 1. For an iV-body system, the implied summation 
would be over i,j — l,3iV. Moreover, if the interactions pro- 
ceed through the usual Newtonian law with no direct impacts 
\7 2 W — 0. If we now label by a, b and c the particle numbers 
(which run from 1 to N) and by k and I the three Cartesian 
coordinates of a particle, it is straightforward to obtain the 
following expressions for the derivatives of W 



Wi 



_dw _8W _ \- r K ac 
dxi drl 2-^ rl, ' 



Wa = 



d 2 W 



d 2 W 



dxidxj drfdr'a 



if a 7^ b and 



d 2 W 

dxidxj 



d 2 W 
dr£dr a 



Ski 



3r nr . 



if a = b. In these equations 



2 

Tab 



I 1 -.2 . / 2 -.2 . / 3 n2 

(r ab ) + [r a b) + (r a b) 



(24) 



(25) 



(26) 



(27) 



and r k ah = r c 

The practical procedure of implementing the criterion de- 
scribed above will therefore consist of using the position and 
velocities of particles for the configurations of the system un- 
der study to obtain u, W, and the quantities defined in (25). 
These are substituted into (23) to find r u . In this way, one can 
calculate the Ricci curvature for various regions of the config- 
uration manifolds of different systems. If the Ricci curvature 
is found to be predominantly negative, we then use the above 
definition to classify systems according to their local stability 
properties. 

There are many advantages to the above setup. For exam- 
ple, what is studied here is the normal deviation due random 
perturbations of trajectories with the same geodesic velocities 
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|| u || = 1 on M and not the deviation of temporal states in 
the direction of maximal growth. Therefore what can be termed 
the 'chaotic pendulum problem' is avoided. In fact, it can be 
seen from (23) that all systems possessing only one degree of 
freedom (those with 3JV = 1) have an r u — at all times. In 
the terminology of classical stability theory it is said that the 
negativity of the Ricci curvature measures the orbital stability 
as opposed to the more strict Liapunov stability (Pars 1965). 
Because of these properties the Ricci curvature is unlikely to be 
negative for multidimensional integrable systems in virial equi- 
librium. This assertion has been checked for the special case 
of the two body problem with circular orbits (where the virial 
relation is satisfied) but obviously needs further investigation. 
For a system in a full statistical quasi-steady state a negative 
Ricci curvature must mean that the system is chaotic and mix- 
ing for all initial conditions since r u is constant for systems in 
a steady state. However, in general, the negativity of the Ricci 
curvature on most of a system's trajectory does not guaranty 
that it is chaotic but only that there is a probability of this 
being the case (the probability increasing with the fraction of 
time spent in the negative region). 

Other advantages of this method are that the integration 
of a large number of linearized equations is avoided and space 
averages on compact manifolds can be compared to time aver- 
ages (e.g., Casetti & Pettini 1993), thus allowing one to check 
for properties such as ergodicity for single particle orbits and 
specially constructed JV-body systems for which M is compact 
(such as those consisting of softened particles enclosed enclosed 
in boxes: e.g., Lynden-Bell 1972). 

3. Applications 

3.1. Some applications that take advantage of the geometric 
setting 

There are at least three applications that can make use of the 
geometric method described above. 

1. Studying the instability properties of individual orbits in 
fixed potentials with compact phase-spaces (in this case 
the formulas in (25) will, of course, have to be modified 
accordingly). 

2. Since r u in Eq. (23) depends only on functions that are 
given by sums over the particle positions and velocities, it 
is constant (up to l/y/N fluctuations) for systems in sta- 
tistical equilibrium. In this case, it is therefore possible to 
replace time averages with phase-space averages (that is 
different realization of same density and velocity fields). 
This time however the averages are made over the full 6iV 
phase-space. The Ricci curvature therefore would give us a 
powerful tool of exploring that space at and around equilib- 
rium solutions. Important questions such as the degree of 
chaos in a system, mixing time-scales, and the variation of 
these properties with particle numbers can then be tackled 
in the 6N phase-space without making assumptions about 
the particle particle correlations. It is very important to 
compare such predictions with those of single particle inte- 
grations in fixed potentials to evaluate the role of discrete- 
ness in the evolution of gravitational systems. For as we 
shall see in the next section, systems with large softening 
parameters appear to have very different 6N phase-space 
structures compared to those composed of point particles. 

3. One can apply the method directly to the results of TV-body 
simulations. The formula for calculating the Ricci curva- 



ture should be easily incorporated into large- N codes such 
as the TREECODE. This would help in interpreting the 
results of iV-body simulations and lead to classification of 
galaxies according to their dynamical instability properties. 
Moreover, one can then study the direction of evolution of 
instability properties of realistic gravitational systems. For 
example it has been argued by Gerhard (1985) that ellip- 
tical galaxies start from chaotic states and evolve towards 
progressively more regular states which are then for long 
times indistinguishable from those of integrable systems. 
An alternative scenario is closer to the conventional dynam- 
ical interpretation of statistical mechanics. In this picture, 
a quasi-steady state is achieved when a system, although 
highly mixing, keeps its macroscopic parameters constant 
and is stable against perturbations to its statistical prop- 
erties. The point being that, if the instability is present for 
most initial conditions of interest (that is if we have con- 
stant negative r u for long enough times), that would mean 
that the system is free to move in the region of interest and 
will tend towards a more probable state. This state would 
(by definition) contain more microscopic states compati- 
ble with it and the system would be free to move between 
them. This situation is more compatible with the interpre- 
tation of violent relaxation as leading to most probable end 
states compatible with a set of constraints (Saslaw 1985) 
since regular states cannot be very probable because they 
lie on 3iV subspaces of the 6N — c dimensional subset of the 
phase-space where all possible states live (c = 10 stands for 
the number of classical integrals used in the reduction of 
the iV-body problem, e.g., Whittaker 1937). 

3.2. Specific application and model parameters 

As a first test we apply the method described above to small N- 
body systems of 231 point particles. The small number enables 
us to integrate the equations of motion with high precision, on 
the other hand the number should be sufficient for the results to 
have some statistical significance. It is essential to check if the 
Ricci curvature method makes adequate predictions about the 
evolution of gravitational systems under these controlled condi- 
tions before applying it to realistic galaxy models where a host 
of auxiliary problems will arise. The small numbers however 
makes it harder to distinguish clearly between the predictions 
of the Ricci method and those of two body relaxation theory. 
A detailed comparison (including dependence of the results on 
N) is better left to another study. 

We choose initial conditions in which the particles are ar- 
rayed into two sheets. An upper one with 11 lines consisting of 
11 particles each and a lower one composed of 11 lines with 10 
particles each. The lines of the upper and lower sheets are po- 
sitioned in such a way that a line in the lower sheet lies at half 
the distance (in the plane of the sheet) between two lines in 
the upper sheet, so as to avoid direct contact of particles from 
two different sheets when the system evolves. The separation 
between the two sheets is taken to be equal to half the separa- 
tion between lines in the same sheet (see the t = snapshots 
in Fig. 3 for example). Such configurations are artificial and 
are therefore likely to quickly and visibly evolve, hence saving 
us the trouble of long time integrations. 

We use units in which the gravitational constant is unity, 
the masses of all particles are also taken as unity. Time in 
these units will be referred to as "physical time" . Three scales 



7 



are used to determine the separation between adjacent par- 
ticles d — 100 x Scale, where Scale takes either the value 
of 1, 10.8 or 100. In what we may call the "main models" 
we give the initial configurations a rigid body rotation corre- 
sponding to an angular momentum of 44275 units around the 
Z-axis (with zero initial velocities in the Z direction), or a ran- 
dom number generator is used to fix the X-Y velocities which 
give rise to a small angular momentum of 62.61 units. The 
three cases of Scale = (1, 10.8, 100) correspond to energies of 
(—41.471, —5.705, —0.635) with corresponding initial virial ra- 
tios of (0.695, 0.0644, 0.00695) respectively. Alternatively, some 
runs are started with a fixed initial virial ratio of 1 and with 
the same energies as above. The angular momentum is adjusted 
accordingly. We shall call these the "equilibrium models" al- 
though they do not start from a detailed dynamical equilib- 
rium. 

The integrations are performed using a variable order vari- 
able step size Adams method as implemented in the NAG rou- 
tine D02CBF using a tolerance of 10~ 13 . The energy is ac- 
cordingly conserved to better than 10 digits (usually 12) for 
a few dynamical times. This accuracy is necessary for a first 
numerical test to eliminate the factor of serious numerical er- 
ror from the interpretations of the results. For the same reason 
we integrate the equations for a relatively short time of about 
four and a half dynamical times (td =(mean density)" 1 ^ 2 ), 
4.24td ~ r b ) corresponding to about 2200 time units for the 
scale = 1 case (in a few cases we have performed longer time 
integrations). The un-softened force law is used, and softening 
is only introduced to study its effect. 

4. Results of numerical experiments 

4-1. Averaging 

In order that our criterion may be useful, the Ricci curva- 
ture must remain mostly negative throughout the evolution 
of a given system. This condition cannot be generally satis- 
fied because r u , as given by (23), is singular as r ab — > for 
any a and b. The time for which r a b is near zero must also be 
small compared to the dynamical time-scale of the system (a 
bounded circular orbit of two particles of separation 0.1 has 
a period of ~ 2 x 10~ 4 td with td the dynamical time of the 
Scale — 1 main model). Therefore, during an encounter of two 
or more point masses, their contribution will fluctuate violently 
and dominate r u , this contribution will have a positive average 
since for a time it mimics that of a two body system. Thus we 
expect r u to be a very "bumpy" function of time for iV-body 
systems of point particles. This is indeed the case, as we can 
see from Fig. 1A, where the Ricci curvature is plotted as a 
function of dynamical time for the rotating main model with 
Scale — 1 (see section 3.2). We will therefore have to try to 
eliminate this bumpiness in the hope that the residual curva- 
ture is negative if the system is mixing. The approach we will 
use here hinges on a theorem by Bogoliubov (P93 and the refer- 
ences therein) which states that for equations of the type (22) 
averaged over small time-scales, the solution is similar to that 
of the original equation averaged over that time-scale if one ex- 
cludes the possibility of parametric instability (e.g., ARN) and 
direct singularities; this is the method we try in this study. We 
will not attempt here to find out what the theoretically opti- 
mum way of doing that should be, but will proceed empirically. 
As a first approach we will just average over a time-scale con- 



taining many "bumps" but still small compared to the total 
integration time. 

Fig. IB shows the time series in Fig. 1A averaged over steps 
of (1/100) x AAtd (corresponding to 22 data points in Fig 1A) 
we see that now the series becomes much more regular. To 
make further progress we remove the contributions from the 
most extreme peaks which can dominate the average giving 
misleading results (these peaks represent points nearer to sin- 
gularities). It is clear from Fig. 1A that peaks with absolute 
values larger than 2 — 4x 10~ 4 are rare and isolated, one there- 
fore is tempted to filter the results so as not to include in the 
calculation of the average any peaks larger than a threshold 
of this order of magnitude. Fig. 1C shows the behaviour of r u 
when a threshold of 2 x 10~ 4 is taken. The striking result is 
the apparent regularity and negativity of the Ricci curvature 
in this case. 

It is now important to check that these results do not sen- 
sitively depend on the values of the threshold used or the aver- 
aging interval taken. Fig. ID shows the behaviour of r u when a 
threshold of 4 x 10~ 4 and 10xl0~ 4 (instead of 2xl0~ 4 ) is used. 
The results clearly appear to be qualitatively similar. Indeed, 
it has been found that the results start to become radically 
different only for a threshold ~ 20 x 10~ 4 . Unless stated other- 
wise, it is implied that it has been checked that the results do 
not sensitively depend on the threshold taken and that a value 
of 2 x 10~ 4 has been adopted. It has also been checked that the 
results do not sensitively depend on the size of the subdivision 
intervals (as long as they are not too small of course). We will 
take averages over a hundred subdivisions of the integration 
interval unless otherwise stated. 

4-2. Rotating versus Non-rotating models 

Rotation, plays a central role in stellar dynamics and is one 
of the parameters that vary along the Hubble sequence. Intu- 
itively, systems where random motion dominates are expected 
to mix faster (at least initially) than ones where ordered ro- 
tational motion does. Observations of elliptical galaxies also 
show them to be highly mixed and relatively (in comparison 
with spirals) relaxed objects so that an effective mixing mech- 
anism must have been at work at some stage (and could still 
be). Our first application will therefore be the comparison of a 
rapidly rotating system to one without (almost) any rotation. 

In Fig. 2 we have plotted the time series of r u for the 
two main models with Scale — 1. Shown by the solid line is 
the series corresponding to the random initial velocities, while 
the dotted line represents the evolution of r u for the rotating 
case. Initially, the random system is much more unstable as ex- 
pected, while later on it starts evolving on a longer time-scale 
(characterized by the small r u ) while the reverse is true for the 
rotating system. 

Explanation of this behaviour of r u may be obtained by 
looking at the spatial evolution of the two systems. The random 
system evolves much faster initially and quickly loses memory 
of its initial configuration in the x — y projection while becom- 
ing diffuse for times beyond t = 1.2 — 1.6 and starts to puff 
up in the z-direction (Fig. 3) thus evolving to a more isotropic 
state. However since it is becoming more diffuse it now takes 
longer to evolve because its density is decreasing. 

On the other hand, the rotating system keeps its shape al- 
most perfectly intact until t = 2 when a distinctive 4 armed 
clustering pattern (Fig. 4) starts to appear accompanied by 





Fig. 1. Ricci curvature time series for rotating "main model" Scale = 1 system: A Unaveraged time series, B with averages taken over a 
hundred intervals, C as in B but with a "filter" of 2 X ltT 4 , D as in B but with a filter of 4 X 1(T 4 and 10 X 1CT 4 (dotted line) 




Fig. 2. Ricci curvature time series for main model with random 
initial x-y velocities and scale = 1 (solid line) averaged as in the 
plot in Fig. 1 C (reproduced here by the dotted line) 



the formation of high density areas. This explains the large 
negative values of the Ricci curvature which in this case rep- 
resents the collective instability. Since geodesic instability, as 
stressed earlier, is independent of the particular evolutionary 
phenomenon in question, this type of relaxation is also included 
in the description. When one interprets relaxational phenom- 
ena in terms of thermodynamic "more probable states" there- 



fore, it should be clear that those types of instabilities could 
also occur. If they do not lead to complete detachment of a 
system into components however, the picture that emerges, is 
that while these phenomena may be important for systems for 
which the classical theory holds, they cannot characterize the 
long term evolution of a generic system and that these inhomo- 
geneities themselves produce chaos and evolution towards more 
probable states under small dissipative or conservative pertur- 
bations (Pfenniger & Norman 1990; Hasan et al. 1993; Pfen- 
niger & Friedli 1991; Friedli & Benz 1993; Zhang 1996). They 
are therefore not typical properties of gravitational systems 
but transient ones. That, roughly speaking, means that regions 
of phase-space where this type of instability occurs should be 
small compared with those where the instability leading to 
more isotropic states arises. 

We now calculate a rough evolutionary time-scale corre- 
sponding to the plots in Fig 2. We take it to be the average 
exponentiation rate for small random perturbations normal to 
the phase-space path of the N-body systems. From the consid- 
erations of Section 2.2 this is given by 



V -2f u J W 



(28) 



where a bar denotes a time average. Now, r u averages to about 
— 1CP 5 over the whole interval (for both systems) while W 
starts at a value of 22 but then rises to about 42, therefore 
we take a value of 30 as a rough mean. This gives a mean 
exponentiation time-scale of about 200 units — or 0.4ro over 
this period of time (0 to 4.4td). 



Fig. 3. Spatial evolution of the model corresponding to solid line 
plot in Fig. 2 

Now, what does this time-scale mean? To answer this ques- 
tion it is important to remember that while the exponential so- 
lutions are local and arise from a linearization of the dynamics, 
this instability exists everywhere because of the negativity of 
the average of the Ricci curvature. Moreover, the instability is 
that of the flow in the full 67V phase-space and not just parti- 
cle orbits. Under these conditions, the divergence will probably 
lead to mixing and filamentation which will radically alter the 
local structure of the 67V phase-space (as mapped by an ar- 
bitrary initial partition of that space). Up to resolution e, the 
phase-space will be completely modified after a time appar- 
ently given by (e.g., GS): 

t £ = lne -1 -r e . (29) 

For a numerical application in double precision the smallest 
resolution available is e ~ 1CP 14 (incidentally, this number is 



Fig. 4. Spatial evolution of the model corresponding to the dotted 
line plot in Fig 3 

not too different from the ratio of the size a star to that of a 
galaxy) so that after a time 

Te ~ 32r e (30) 

there is complete change in the system's phase-space for all 
practical purposes. For an ensemble of trajectories having the 
same energy and integrated independently in a fixed potential, 
Merritt & Valluri (1996) found that complete mixing did in- 
deed occur over a number of exponentiation time-scales (mea- 
sured by Liapunov exponents in their case) comparable to that 
given by (30). 

Examples of systems that were proved to behave in this 
manner include the general class of C-systems (e.g., Arnold's 
famous cat map: LL). Although we cannot claim that the pro- 
cess is as rigorous in our case, with its non-compact phase-space 
and sign indefiniteness of the 2d curvatures, we will follow the 
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assumption that it is at least qualitatively similar (the negativ- 
ity of the averaged Ricci curvature, the results of GS and Kan- 
drup 1990a,1990b suggest that this may be so). In that case, we 
have complete modifications of our systems on a time-scale of 
about 13 dynamical times. Obviously, our two examples have 
evolved significantly over a time-scale of this order of mag- 
nitude — albeit each in different ways and non-uniformly in 
time. We shall see in Section 4.4 that on a time-scale of 12td 
there is indeed complete modification for one similar system. 

Finally we note that formula (2) predicts an exponential 
time-scale of about I.Itd which is somewhat larger than the 
one derived from the Ricci curvature but within reasonable 
bounds considering it is only an order of magnitude estimate 
derived for infinite systems and does not take into account 
effects due to the large scale gravitational field (see Kandrup 
1989,1990b) that lead to the clustering in the rotating case for 
example. 

4-3. Varying the energy and virial ratio 

We now consider systems with lower densities while keeping the 
angular momefinger ntum constant. For systems initially in a 
state of collapse (vir < 1) this amounts to varying the energy 
— making it larger. For these systems, the initial virial ratios 
(and hence the velocities) is very small, therefore the results 
for the rotating systems are more or less similar to those for 
systems starting with random X-Y velocities. We will therefore 
focus on the latter type only. 














Fig. 5. Ricci curvature time series for random main model systems 
with: A Scale = 10.8, B Scale = 100 averaged as the plot in Fig. 1C 
(reproduced here by the dotted lines) 

In Fig. 5 we compare the Ricci time series for initial condi- 
tions with Scale = 10.8 (Fig. 5A) and Scale = 100 (Fig. 5B) 
(shown by the solid lines) with those of the Scale = 1 (dot- 
ted line) . Clearly, there is a trend towards larger time-scales of 
evolution for the less dense systems, especially if one takes into 
account that the ratio of the kinetic energy W at virial equilib- 
rium of the Scale = 1 system to that of the Scale = 10.8 and 
the Scale = 100 systems is 7.3 and 65.9 respectively. These 



results may seem counterintuitive since one knows that for ini- 
tial conditions far from dynamical equilibrium there is a well 
known instability at work — namely violent relaxation. The 
confusion is resolved however when one notices that the dy- 
namical time for these two systems is much larger than the 
original — being 35.49 and 1000 times larger. Thus, based on 
the above considerations, we expect a system characterized by 
the Ricci curvature shown by the solid line in Fig. 5B to evolve 
extremely slowly in "physical" time units but relatively fast, 
compared to the one shown by the dotted line, in terms of dy- 
namical time. The fact that the kinetic energy varies rapidly 
along the motion prevents us however from calculating any 
precise estimates. 



10 5 HJ 5 



10 s HO 5 



Fig. 6. Spatial evolution of the model corresponding to the solid 
line plot in Fig. 5 B 

By looking at Fig. 6 one can clearly see that this is the 
case. With every interval between "snapshots" in this figure 
corresponding to a thousand intervals in Fig. 3 one sees that the 
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Scale = 100 system evolves much more slowly in physical time. 
In terms of intrinsic dynamical times however, it is clear that it 
evolves much faster than the Scale = 1 system; having quickly 
lost all trace of its initial spatial configuration it then appears 
to be evolving towards a state characterized by an increasingly 
tight core and diffuse "halo". At this stage one notices, the 
almost constant (drifting slowly towards lower) negative values 
of r u in fig 5B (solid line) beyond t ~ 3, showing that the 
system is becoming more and more mixing, in agreement with 
discussion in section 3.1. 

4-4- Systems starting in virial equilibrium 

We now keep the energies constant and increase the angular 
momenta of the rotating systems in such a way as to have 
virial equilibrium in the initial state. This will be accompa- 
nied by contraction in the scale of the system. For example, 
the Scale — 100 equilibrium models will have an inter line 
spacing of about half that of the main models with the same 
energy. The systems with random x — y velocities start from 
the same spatial configurations as the rotating ones and with 
the absolute values of the velocities rescaled accordingly. 
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Fig. 7. Ricci curvature time series for random main model systems 
with: A Scale = 10.8, B Scale = 100 averaged as the plot in Fig. 1C 
(reproduced here by the dotted lines) Ricci curvature time series of 
the rotating equilibrium models (solid lines) compared with the main 
(collapsing) models: A Scale = 1, B Scale = 10.8, C Scale = 100. 
In B and C averages are taken over 50 time intervals (instead of a 
100) 



For a system in a steady state we expect r u to be constant 
in time. If, in addition, the motion is regular (e.g., our sheets 
rotate rigidly) we expect r u to be very near zero. Nevertheless, 
due to the differential rotation which is accompanied by collec- 
tive instability and clustering, this state of affairs cannot last 
too long. Similarly, due to the instability of the initial state 
in the random case, the early equilibrium is lost. These con- 
siderations are confirmed by plots of the Ricci curvature (Fig. 
7) where values for the equilibrium cases (solid lines) are com- 
pared with those of main models of the same energies (dotted 
lines). Although the curvature is still negative even at early 
times, the evolution time-scale given by Eq. (30) is now about 
70 dynamical times. It also is immediately apparent that the 
more distant the original system was from virial equilibrium 
the bigger the difference of its initial behaviour from its equi- 
librium counterpart. For example, for the Scale — 1 system 
(Fig 7A) the dotted and solid lines are almost indistinguish- 
able, while for the other two systems the contrast in the early 
evolution is much clearer. Later on however, non-equilibrium 
systems quickly tend towards virial equilibrium, making the 
contrast weaker. (The reason for the positive bumps in the 
time evolution of the curvature for these models is that, de- 
spite the fact that we average over 50 time intervals, the ir- 
regularities cannot be suppressed because the filtering is too 
large the fluctuations now are much smaller because of the 
lower density.) It is clear from these graphs that, because the 
dynamical time of the equilibrium systems is smaller than the 
collapsing ones (about a third for the Scale = 100 models), the 
equilibrium models actually evolve slower than the collapsing 
systems in terms of intrinsic dynamical times throughout most 
of the evolution. From that, and from the results of the previ- 
ous subsection, we conclude that systems starting from a state 
of collapse are possibly more unstable than ones starting near 
virial equilibrium. The difference does not however appear to 
be very large. 

The x — y spatial evolution in time of the rotating Scale = 
100 equilibrium model is shown in Fig. 8 where the integration 
was carried over a time interval corresponding to the physical 
time of integration of its main model counterpart with same 
energy. It is clear that by t = 4.48 the clustering is at least 
as pronounced as in Fig. 4 at t = 4.4 (a similar comparison 
can be made with the corresponding equilibrium model which 
behaves in a very similar manner). This is not surprising given 
the scale free nature of gravitational interactions. It is therefore 
interesting to see what the rough time-scales derived from the 
Ricci curvature might tell us about this apparent coincidence. 
Let us label the quantities relating to the Scale = 1 equilibrium 
initial state by indices 1 and those relating to the Scale = 100 
equilibrium system by 100 . The ratios of the two evolutionary 
time-scale are then 

ti y r ul00 WlOO TD100 

Both systems are always near virial equilibrium throughout the 
evolution therefore the W's are constrained to be near values 
given by W = -E so that Wi ~ 41.53 and W wo ~ 0.63, the 
averages f u i and f u ino are calculated over the interval t = 
to t = 4.4 to be f u i = -1.6 x 10~ 5 and f u ioo = -5.8 x 10~ 7 . 
Finally 

TDioo = / droo \ 3/2 = (5Q 17) 3/2 = 355 41j (32) 

TBI \ Ol / 
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Fig. 8. Longer time behaviour of the Ricci curvature and the cor- 
responding x-y projection of the spatial evolution for the system 
corresponding to the solid line plot in Fig. 7C. Filtering is taken 
equal to 2 X 10 -5 and averaging is done over 150 intervals. 



where di and dioo are the initial separations of adjacent parti- 
cles on the sheets. Inserting these numbers we get 



about 12 dynamical times as expected from Eq. (30). To il- 
lustrate how the results depend on filtering, we have chosen a 
filtering here of 2 x 1CP 5 instead of the value 2 x 1CP 4 usually 
used. Hence the smoothness of the time series up to t — 4.4 in 
comparison with the one in Fig. 7C. The positive bump on the 
right is due to the filtering becoming too small relative to the 
absolute value of r u so that the value in this last interval has 
little statistical significance. 

Similar results are obtained for the cases with random ini- 
tial x — y velocities. Here, the formula for two body relax- 
ation for a Maxwellian distribution (BT formula 8.71 which 
is more or less appropriate in this case) gives a ratio of 
tioo/t1 = 0.87 3 = 0.66 for these two cases. 

4-5. The effect of softening 
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: 0.97. 



(33) 



Now, considering that the systems differ radically in energy, 
angular momentum and dynamical time-scales this result is 
rather impressive — especially if one recalls the rough way in 
which we have averaged and filtered the data. It may be useful 
to note that Eq. (2) predicts the same behaviour, we have 



W 1/2 W 1/2 



,1/6 



t d ~ {Scale x W) 1/2 



which gives 
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(34) 



(35) 



Looking at the longer time behaviour in Fig. 8 we see 
that the clustering continues leading to almost total separa- 
tion of the components. The corresponding behaviour of r u 
is also shown. As would be expected, the clustering gives rise 
to shorter evolutionary time-scales. Three characteristic time- 
scales can be isolated: from t — to t — 7 the system can still 
be considered to be one part and is connected until t ~ 10 but 
then detaches after that. These stages are clear from Fig. 8 . 
Thus there is complete modification of the initial state after 




Fig. 9. The effect of softening: Same as in Fig. 8 but when the force 
law is softened: A When the softenening radius is equal to 10% the 
inter-particle distance on the same line, B same as in A but when 
the softening radius is equal to 80% the initial intcr-particlc distance 
on the same line. The x — y coordinate evolution (top) correspond 
to the system in B 
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Fig 9A shows the un-averaged Ricci time series for a sys- 
tem starting with the same initial conditions as the one in 
Fig. 8 but when the potential contains a softening parame- 
ter (BT formula 2.194) of 10% the initial separation between 
particles in the same line, while Fig. 9B shows the the evolu- 
tion of r u for the same system but when the softening radius 
is about 80% of the original separation between particles in 
the same line. First, one notices the absence of large fluctua- 
tions in the value of r u . This is to be expected since formula 
(23) no longer contains singularities at close encounters. The 
second and more fundamental effect however is that the Ricci 
curvature is now positive. This is due to the fact that the last 
term on the right hand side is now non-zero and is always pos- 
itive since it corresponds to the density. For small softening 
(say 1% of inter-particle distance) the original behaviour of 
r u is recovered, however as the softening increases this term 
becomes large. In addition, the second term in (23), which in 
general gives a highly fluctuating contribution the average of 
which is small, starts becoming positive too. 

The above means that in the presence of significant soft- 
ening, N-body gravitational systems no longer approximate the 
uniformly hyperbolic Anosov C-systems because the Ricci cur- 
vature is no longer negative on average. That is, there is a 
fundamental change of structure in the phase-space of iV-body 
systems when softening is introduced. Therefore if one takes 
the view that a softened iV-body system is an approximation 
of the continuum formulation of a problem, one must realize 
that this formulation has solutions which can behave in a dif- 
ferent manner from those of a discrete system. Indeed it was 
found that the evolution of the softened systems differed con- 
siderably from the un-softened ones. For example in Fig. 9 we 
also show the x — y evolution from the same initial conditions 
as in Fig. 8 but when a softening parameter of 80% the original 
separation of particles is present in the force law. Instead of the 
clustering pattern that was obtained in the un-softened case, 
there is now a type of filamentary structure taking its place 
and the system never completely detaches into components. 
(In fact when the simulation is continued these filaments are 
destroyed and one gets a nearly isotropic system.) Although 
this system is still unstable, the dominant mechanism of insta- 
bility is likely to be different. It was observed for example that 
while in the un-softened system the virial ratio stays close to 
1, in the softened case it departs significantly from that value, 
being less than 1 for the time-scale of the plot then increasing 
significantly later. This may suggest that the mechanism of in- 
stability here is a parametric instability related to that effect 
(Cerruti-Sola & Pettini 1995). For systems with intermediate 
softening a mixture of different mechanisms could be present. 
When the negativity of the curvature is dominant, a system 
approximates a discrete one. This criteria may be useful in 
choosing softening parameters (or laws) in numerical simula- 
tions that eliminate the effect of close encounters but do not 
qualitatively change the results. 

The difficulty we are encountering here however is a fun- 
damental one. As long as the system we are considering is dis- 
crete and direct collisions are negligible the offending term in 
Eq. (23) is zero and the Ricci curvature is likely to be negative. 
As mentioned in the introduction and in the beginning of sec- 
tion 2.3, this is even more likely to be the case as N increases. 
On the other hand the smoothing of the density distribution 
leads to positive curvature. The structure of the phase-space 
(which is the cotangent bundle of the configuration manifold 



the curvature of which we are calculating) in the two cases is 
then presumably radically different. This suggests that it is 
not obvious that large-N but discrete systems are well approx- 
imated by the continuum limit even if the discreteness noise in 
the force decreases as N increases. There appears to be a type 
of "phase transition" in the structure of the 6N phase-space 
when passing from one limit (large- N but discrete system) to 
the other (continuous system). This is the same conclusion we 
were led to from the general considerations of the introduction. 

Obviously these issues will require much more thorough 
investigation which is beyond the scope of this exploratory pa- 
per. We will however briefly comment on this problem from 
the perspective of the scalar curvature which we now discuss. 

4-6. Scalar curvature 

When the two dimensional curvatures are averaged over all 
geodesies that originate at a point of M as well as all possible 
directions normal to the geodesies, one obtains the curvature 
scalar 

i? = ^fc u , n = ^r u . (36) 

u.n u 

With the quantity R/3N(3N - 1) replacing r u /(3N - 1) in 
Equation (22) one can define a corresponding instability cri- 
terion and an associated time-scale. This is obviously a more 
drastic approximation, and because R will not depend explic- 
itly on the velocities (since we have summed over the directions 
in the tangent space TM) it cannot contain generic informa- 
tion about a system and its use for short time characterization 
of the dynamics can be especially dangerous (e.g., Cipriani & 
Pucacco 1994). Nevertheless, the similarity of the third term 
in the formula (23) for the Ricci curvature to the expression 
for R (GS formula (27)) suggests that when this term is dom- 
inant, the scalar curvature may give a reasonable description 
of the dynamics. The fact that estimates of evolutionary time- 
scales resulting from the order of magnitude formula of GS 
agreed within reasonable bounds with the ones obtained from 
the computations of r u suggests that this term may be impor- 
tant in some cases since that formula is based on an estimate 
of the scalar curvature. 

To see if this is true, we have computed the scalar curva- 
tures along the motion of the systems described in the previous 
sections. We have found that using a hundred averaging inter- 
vals and filtering equal to twice the absolute value of the max- 
imum scale in the figures gives reasonably smooth time series 
for R that do not sensitively depend on the averaging and the 
filtering threshold used. Fig. 10 shows the resultant time series 
for the rotating equilibrium models. Clearly, those results agree 
qualitatively with those inferred from the behaviour of the cor- 
responding time series for r u in Fig. 7 (although differing in 
some detail). The time-scales also agree (when R is divided by 
3N(3N — 1) as described above). This agreement is because 
out of the three non-zero terms in Eq. (23), the first is small 
unless the velocity vectors of particles are almost aligned with 
their acceleration vectors and the second is highly fluctuating 
and averages to very small values (this is because it describes 
the variations of the force at the particle positions which, as 
one would expect, depends sensitively on the position of close 
neighbours). This leaves the third term which is the one that 
appears in the formula for the scalar curvature. 
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Fig. 10. Scalar curvature for the rotating equilibrium models A 
Scale = 1, B Scale = 10.8, C Scale = 100 



The scalar curvature — like the Ricci curvature — also 
takes positive values when significant softening is used. As 
was noted in the introduction, the continuum approximation 
with a steady state potential reduces the TV-body problem to 
N one-particle problems in a given potential. In this context, 
the above discrepancy between the value of the curvatures of 
softened and un-softened cases may perhaps be understood by 
noticing that the scalar curvature for N = 1 is always positive 
no matter what the potential is (assuming positive density). 
The 6-dimensional phase-space is therefore far from approxi- 
mating a hyperbolic structure like that of C-systems. Even in 
this case, chaos may still occur (some two dimensional curva- 
tures may still be negative and there is also the possibility of 
parametric instability: ARN; P93) but it will not likely be as 
robust or widespread as in systems with negative curvature. 

5. Conclusions 

The assumptions of classical non-evolutionary galactic dynam- 
ics may not be satisfied in the presence of significant amount 
of chaos (section 1). This makes it important to develop and 
to test methods characterizing such chaotic behaviour. Geo- 
metric methods have the double advantage of being local and 
comparing normal deviations between trajectories of dynami- 
cal systems and not the total divergence of temporal states. For 
higher dimensional systems these properties may be important 



in distinguishing between true phase space mixing which can 
be accompanied by changes in the physical characteristics of a 
system and phase mixing which conserves the action variables. 
Geometric methods are therefore better suited for studies of 
the short time evolution of galaxies than other (more tradi- 
tional) measures of chaos (section 2 and 3.1). 

In large dimensional systems it may not be practical to 
determine the stability of a trajectory to all possible perturba- 
tions. The next best thing is to determine the average diver- 
gence of trajectories due to random perturbations. One way 
of doing this is by examining the Ricci curvature of the La- 
grangian configuration manifold of a dynamical system (section 
2.3). To check the effectiveness of such an approach for iV-body 
gravitational systems we have calculated the Ricci curvature 
for several small iV-body systems integrated with high preci- 
sion (section 3.2 and 4). The results of these experiments show 
that: 

1. When properly averaged to get rid of the contributions of 
close encounters the Ricci curvature is almost always neg- 
ative, confirming that gravitational systems are unstable 
(e.g., Miller 1964; Goodman et al. 1993; Kandrup et al. 
1994) and that the main mechanism of instability is the 
negativity of the curvature of the configuration manifold as 
predicted by Gurzadyan & Savvidy (1984,1986) and Kan- 
drup (1990a,1990b). 

2. The Ricci curvature is more negative (hence predicts 
shorter evolutionary time-scales) when a system develops 
pronounced macroscopic instabilities (e.g., plasma type col- 
lective instabilities). In this case the rates of spatial macro- 
scopic evolution of the different systems both relative to 
each other and in terms of evolution time-scales was rea- 
sonably well described by the time-scales derived on the 
basis of the Ricci curvature calculations. 

3. When expressed in terms of dynamical times, evolutionary 
time-scales appeared to be slightly longer for systems start- 
ing from virial equilibrium than those starting from a virial 
ratio less than one. However, since when the kinetic energy 
varies significantly, chaotic behaviour cannot be described 
by the negativity of the Ricci curvature alone (Cerruti-Sola 
& Pettini 1995) this conclusion remains to be confirmed. 
(However in the case of large iV-body systems near virial 
equilibrium the objections to the use of the Ricci curvature 
as outlined in the aforementioned paper are not likely to 
be important since the second and third term of their Eq. 
(26) are then very small. By using the Ricci curvature and 
eliminating large fluctuations due to close encounters we 
have therefore implicitly assumed that it is the instability 
in that (large- N) limit that interests us and not effects due 
to small scale fluctuations.) 

4. In the presence of significant (but not very large) soften- 
ing, the Ricci curvature becomes positive. This probably 
means that the phase-space structure of softened systems 
is radically different from that of point particles. This has 
consequences for the interpretation of results of numerical 
simulations. More fundamentally however this result may 
be interpreted to mean that there is no continuous tran- 
sition from large-N discrete N-body systems to continuous 
ones. It may therefore explain why large iV-body spherical 
systems are found to approximate exponentially unstable 
C-systems while it is known that motion in smooth spheri- 
cal potentials is separable. More work however is needed to 
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fully understand the meaning of this effect and its possible 
consequences. 

5. Results derived on the basis of the scalar curvature agreed, 
in general, with those obtained from the evolution of the 
Ricci curvature. This shows the instability to be quite a 
robust phenomenon. 

We have not looked at properties such as energy relaxation 
here. As was mentioned in the introduction this can have dif- 
ferent time-scales from that of the instability of trajectories. 
To see how this may be the case we consider the change in 
the Hamiltonian H — T + V of a test particle i in an TV-body 
system. Using the Hamiltonian equations we find that along 
the motion this is given by 

dHj _ dHj _ dVj 

dt dt dt 1 ' 

with the interaction potential Vi given by 

V i = J2 V ij ( 38 ) 

3 

being a function of the positions of the remaining N — 1 parti- 
cles which are of course time dependent functions of the initial 
conditions. The change in energy of particle i along its path is 
then given by 

3 3 

which is the sum of the energy changes due to the individual 
interactions and therefore could proceed on time-scales simi- 
lar to that given in (1). This does not of course mean that 
energy relaxation cannot be enhanced for systems consisting 
of particles with different masses or those out of virial equi- 
librium or where collective motion or large scale inhomogene- 
ity or anisotropy occurs. In these cases chaotic behaviour can 
be important for energy relaxation. Indeed there is some ev- 
idence that in some situations two body relaxation estimates 
are inaccurate even for energy relaxation (El-Zant 1996a). In 
general however there may be phase-space "barriers" across 
which diffusion is slow. This may prevent some quantities from 
relaxing even when chaos is present (discussions of the issue 
of phase-space transport in higher dimensional Hamiltonian 
systems can be found in Wiggins 1991 or Benettin 1994). Nev- 
ertheless, chaos implies exponential divergence normal to the 
phase-space trajectory leading to diffusion in at least some of 
the action variables which determine the physical characteris- 
tics of a system (as opposed to phase mixing which conserves 
the action variables). It therefore has important consequences 
for the behaviour of dynamical systems as was stressed in the 
introductory section of this paper. This has long been gener- 
ally recognised in various branches of Physics and Mathematics 
(e.g., Sagdeev et al. 1988) but not always in stellar dynamics. 
The question therefore is how widespread is the chaotic be- 
haviour in realistic N-body realizations of the different galaxy 
types and what are the exponentiation time-scales. This is what 
we hope to find out using the method examined in this paper. 
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